function [rm, w] = rmean(d, dim)
% rmean  - re-weighted mean
%
% FORMAT:       [rm, w] = rmean(d [, dim])
%
% Input fields:
%
%       d           data array
%       dim         optional dim (otherwise first non-singleton)
%
% Output fields
%
%       rm          re-weighted mean
%       w           weight

% Version:  v0.7f
% Build:    8110522
% Date:     Nov-05 2008, 10:36 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~isnumeric(d)
    error( ....
        'BVQXtools:BadArgument', ...
        'Invalid argument supplied.' ...
    );
end
if nargin < 2 || ...
   ~isa(dim, 'double') || ...
    numel(dim) ~= 1 || ...
    isinf(dim) || ...
    isnan(dim) || ...
   ~any(dim == 1:ndims(d))
    dim = findfirst(size(d) ~= 1);
    if isempty(dim)
        dim = 1;
    end
end

% compute general mean first
sz = size(d);
ne = sz(dim);
sz(dim) = 1;
rs = osz(sz);
rs(dim) = ne;
ed = 1 / ne;
rm = ed .* sum(d, dim);
ow = ones(rs);

% now iterate until difference to last iteration gets small
mi = 50;
mb = sqrt(1 / ne);
while true
    w = (d - repmat(rm, rs)) .^ 2;
    for dc = 1:ndims(d)
        if dc ~= dim
            w = sum(w, dc);
        end
    end
    w = 1 ./ sqrt(w);
    w = min((1 / mean(w)) .* w, 2);
    nm = (1 / sum(w)) .* sum(repmat(w, sz) .* d, dim);
    wd = abs(w - ow);
    if mi < 1 || ...
        sum(wd) < mb
        rm = nm;
        break;
    end
    rm = nm;
    mi = mi - 1;
end
